Cavitation induced by explosion in a model of ideal fluid 
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We discuss the problem of an explosion in the cubic-quintic superfluid model, in relation to 
some experimental observations. We show numerically that an explosion in such a model 
might induce a cavitation bubble for large enough energy. This gives a consistent view for 
rebound bubbles in superfluid and we indentify the loss of energy between the successive 
rebounds as radiated waves. We compute self-similar solution of the explosion for the early 
stage, when no bubbles have been nucleated. The solution also gives the wave number of 
the excitations emitted through the shock wave. 

Contents 



introduction. 



II The Gross-Pitaevskii equation, 



rn The modeJI 



iV The explosion 



V Decomposition of the process 



Vi Analysis of the exponents^ 



Vil Renormalisation Cjlroup calculations 



VIII Conclusion- Acknowledgments 



2 
3 
6 
7 
9 
11 
16 
19 



1 



I. INTRODUCTION. 



Cavitation is a physical process involving such aspects as erosion, bubble formation, 
sonoluminescence and first order phase transitions. The motivation for this work comes 
from experiments on superfluid Helium and we will generally speak in this context, 
although the model might be applied to others fluids, and to non linear optics [^. The 
experiments of study the cavitation process in superfluid Helium (He^). There, a 
semi-spherical convergent sound wave is produced in the liquid. At the center of the set- 
up, the superfluid alternates between being compressed and being under tension {i.e. at 
negative pressure). During the tensile strength period, a bubble can be nucleated by thermal 
activation. This bubble acquires a kinetic energy through the negative pressure region so 
that it grows inside the liquid bulk until it reaches a maximal radius (a few hundred /xm 
can be obtained ) determined by its kinetic energy and the mean liquid pressure. Then, the 
bubble collapses under this positive pressure. After its collapse, a secondary bubble, called 
the rebound bubble is observed. The collapse of bubbles is a catastrophic process that has 
been widely studied, particularly because of its ramification in industrials application (for 
a review, see 0). Sonoluminescence arises also during the collapse of the bubble 0. The 
collapse is generally followed by a shock wave that is often a cause of important damages. 
The goal of this paper is to show that in the particular case of superfluid, the rebound 
bubble can actually be nucleated by this shock wave. It has been noticed already that 
detonation in water can nucleate cavitation bubble through the tension shock wave that 
follow an explosion (see 0). 

We will first introduce a model recently used for phase transition in systems like su- 
perfluids Then we will numerically study the problem of explosion in this model as 

the collapse of the bubble might indeed be interpreted as an explosion. The bubble that 
will form beside the shock wave might be understood as the rebound bubble found in the 
experiments. We show also that the explosion process exhibits two distinct regimes: for 
small time, it follows a pure explosion behavior whereas an interaction between sound wave 



2 



propagation and interface dynamics dictates the larger times. We finally examine the large 
energy limit, where we can expand the explosion in terms of a self- similar solution for small 
time. 



II. THE GROSS-PITAEVSKIi EQUATION. 

The Gross-Pitaevskii (GP) equation has been often used as a model of superfiuid [§. 
It describes the time evolution of a complex function '?/'(x;t), called the condensate wave 
function; it reads: 

ihdtij{v, t) = -— V^(r, t) + g\ij(v, t) |2^(r, t) (1) 
im 

where h is the Planck constant, m the mass of the particles (for superfiuid, this is miif.i) 
and g is the strength of the potential. 
One can write: 

and in analogy with quantum mechanics terminology, one can identify p as the particle's 
density and ^0 as the velocity potential: 

v(r,t) = -V0(r,t) 
m 

It follows that the equation is conservative: the total number of particles is conserved by 
the dynamics. 

N = j dv\^\^ 

The dynamics is also hamiltonian which means that one can define an energy H such 
that: 

where ifj* is the complex conjugate of ip and with: 
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The halmitonian structure implies that the dynamics is reversible. 

A set of equations of hydrodynamic form (for p and (p) can be deduced from the GP 
equation: 

dtp + V- (pv) = 

(hct>\ 1 g 

-dt{ — ] = -—^ ■ — ^ + -V +—P 



The first equation is the mass conservation and the second one can be viewed as an 
equivalent to the Bernoulli equation for fluids. Indeed the pressure is divided into two 
terms: the first one, gives a static pressure: P = i^p^ whereas is called the quantum 
pressure term, because it annihilates when /i — > 0. This term is a reminiscence that the 
dynamic is deduced from a Schrodinger equation. Without this term, the equation would 
read exactly as the Euler equation for perfect fluid with a given state equation for P{p). In 
fact, the quantum pressure changes this Euler dynamics, (while conserving the total energy) 
whereas for real fluid, it is the viscosity that stabilizes the flow, with a dissipative dynamics. 

Actually, no damping terms (like viscous term) are present in this dynamics; this is 
necessary for a consistent model of superfluid; also, this equation admits only liquid-like 
solutions (^0 = y^e~*^*, where po is the liquid density, constant); the perturbations 
around such solution: 

respects the dispersion equation: 

2 9- Po 2 , 

m 4m"^ 

cu being the energy of the so-called quasiparticles and p their momentum {p — h ■ k, k 
wavenumber). For low wave numbers, the quasiparticles are phonons, i.e. sound waves with 



sound velocity Cg = w^^. The phonon spectrum for low p has been pointed out to be a 



crucial property needed for modeling superfluids, according to Landau theory |Tl| . For large 
wave numbers, the spectrum corresponds to the one for free particles (the kinetic term is 
dominant). Then, unfortunately the model does not describe the roton part of the superfluid 
spectrum, although a non local potential in equation (|l]) would allow for such a spectrum 
. The crossover between the two regimes occurs for a typical length, called the coherence 
length ^o: 

^0 



Finally, such a system contains by construction another important feature of superflu- 
idity, the so-called quantum vortices. They are in fact topological defects associated to the 
complex order parameter ■?/' and reflect that its phase might be multivalued (modulo 27r). 
As the phase of ip is related to the velocity potential, the circulation of the vortices is a 
multiple of h/m, as predicted by Onsager for superfluid |jl3 . 



Briefly, all these properties have made the Gross-Pitaevskii equation a reasonable model 
for a superfluid at T = Kelvin, and it is often used due to its balance between simplicity and 
sufficient physical ingredients. Therefore it has been particularly used for numerical studies. 
This equation might also be simply considered as a fluid dynamics model, satisfying an 
Euler equation but with an additional term arising from the quantum pressure. This term 
becomes relevant for dynamics on length scales smaller than the coherence length ^q, and in 
particular it stabilizes the vortex core. Such a fluid model presents an alternative view of 
fluid dynamics, in that divergences (shocks...) are eliminated by dispersion rather than by 
dissipation (as the viscosity does for real fluids). Aside from this theoretical point of view, 
such approach might give interesting prospects for real fluids. 

In the context of this paper, it is relevant to introduce a dimensionless form of equation 
(2). It is the so-called non linear Schrodinger equation (NLS), obtained by trivial rescaling 
of the space, the time and the mass; it reads: 

tdt^ = ~A^ + \^\'^ (2) 



Therefore, if po is the mean density of particles, the sound velocity is worth ^/po, the 
coherence length = typical time scale being 1. 



III. THE MODEL 

The equation above describes, in fact, the dynamics of a monophase fluid: precisely 
the model allows only one thermodynamically stable phase, called liquid phase. Therefore 
the model is not relevant for any problems involving liquid-gas transition or any first order 
phase transition. Particularly, all the interaction that can occur in high speed flow between 
vorticity and cavitation are lost in such an approach. It has been shown that a small change 
in the NLS equation can give a consistent first order liquid-gas transition without changing 
the other properties of NLS; it has then be used for studying standard problems like flow 
around an obstacle in two space dimensions and coarsening processes [^lOl- This subcritical 
non linear Schrodinger (SNLS) equation reads in a dimensionless form: 

idtt/j = -^Aij + (2p, - po)po^ - 2p,\ij\^ij + IV-IV (3) 

Pc and Po are respectively the critical and the mean densities as will be explained below. 
The equivalent set of equation for (|) reads: 

dtp = -V ■ (pv) (4) 
5,0 = ^A(P'/') - ^(V0)' -p' + 2p,p - po(2p, - Po) (5) 

They have the same structure than the one deduced for the Gross-Pitaevskii equation. 

The Bernoulli equation allows us to define the static pressure P (forgetting the quantum 
pressure contribution) and therefore a sound velocity might also be derived: 

2 

P = ^p^-PcP^-, (? = Mp-Pc) 

They are shown on figure (||). It appears clearly that pc corresponds to the spinodal 
decomposition density, where the sounds velocity vanishes. As no temperature exists in 
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the model, the pressure dependence in p plays the role of a state equation for the fluid. 
Particularly, with po being the liquid density, one can investigate the pressure difference 
between the liquid phase and the gas phase (which is at zero pressure). For po < fpc 
the liquid pressure is smaller than the gas pressure, which means that the liquid phase is 
metastable relative to the gas one; we have the opposite situation for po < fpc- p = |pc is 
the density for which the liquid and the gas pressures are equal. Therefore, the density of 
the liquid when gas and liquid coexist has to be p = |pc. 

In one dimension, exact solutions of equation are known for a given number of 



particles (see [|lOl)- Consequently, the gas-liquid interface can be exactly computed. It 
connects a region of zero density to a region of |pc density. The energy of such a solution 
gives the surface tension a {i.e. the energy of the liquid/gas interface) of (^: 



a 



16^6 

We have available now a mean field model of a first order transition. The fluid obeys an 
Euler-like equation of motion via a complex order parameter if). Even though the model is 
not entirely physically realistic, it has been used for quantum flows, and we believe that it 
has some important qualitative features. Particularly, as a mean field model, the liquid-gas 
interface is automatically solved by equation (^ without using special treatment. 



IV. THE EXPLOSION 

Using the model introduced above, we will now focus on the particular problem of strong 
explosion. Such questions have been investigated long ago for gases |]14| and interesting self 
similar dynamics have been pointed out. As proposed in the introduction, we will consider 
that the collapse of the initial bubble in the experiment leads to the formation of an explosion 
in the bulk. In SNLS, such an explosion corresponds to a peak of energy (or equivalently 
of density) centered at the origin of the collapse. Numerically, the initial conditions for the 
explosion with a peak of density centered in r = will be taken as: 
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Here, po is the liquid density, which is chosen such that the hquid is the stable equilibrium 
phase (the gas being the metastable one). For now on we will have po = 1 and pc = 0.6 
(therefore all the physical quantities are of order one). 6 is the amplitude of the excess 
pressure due to the explosion and cr is its width. For a strong explosion, where 5^ ^ 1, the 
energy Eq of the explosion is given by: 

Physically, a has to be on the order of few coherence lengths (we will show mainly results 
with a = 1 although we did try a large range of values). Also it appears that the process we 
will describe below is robust and does not depend strongly on the specific initial conditions. 

The numerical simulations have been performed with a finite difference Crank-Nicholson 
scheme [|T5|, that preserves the number of particles exactly at the first order in time. The 
code has been employed for 1, 2 and 3 space dimensions with cylindrical and spherical 
symmetry for 2 and 3 dimensions respectively. We did in fact neglect the Rayleigh- Taylor 
(RT) instability during the simulations. The RT instability appears generally only during 
the collapse and is enhanced by the proximity of a boundary (see |^). However, whether 
the collapse is symmetric or not, it will still give rise to an explosion (perhaps weaker in the 
asymmetric case). 

One of the main numerical limitation comes from the fact that the sound velocity is 
proportional to the local density for large densities p ^ po; therefore, for investigating 
strong explosions, one needs to deal with large sound velocity in the initial density peak. 

Figure (|^) shows the evolution of the density profile in spherical geometry for 5 = 3 
and a = 1. As expected, a gas bubble is nucleated backwards the shock. Also, a train of 
spherical waves is emitted during the process. The bubble grows until a maximal radius is 
reached and then collapse occurs. In figure d), the bubble has just collapsed, giving rise to 
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a secondary explosion, much smaller than the initial one, because part of the energy has 
been transformed in excitations waves. This latest stage gives a convincing proof that the 
collapse of a bubble in such a model might be investigated as an explosion. The bigger the 
initial explosion is, the bigger the maximal radius and the secondary explosion are. Also, 
if the initial energy is lowered, no bubbles are nucleated below a certain critical energy. 
Formally, the secondary explosion can nucleate an new bubble and so on as long as it has 
enough energy, although we did not investigate large enough initial explosions to see at least 
a secondary bubble. 

Whether or not a bubble is nucleated, the explosion in SNLS exhibits a general picture: 
the explosion expands until a maximum radius and then collapse in a secondary explosion, 
which expands again until its maximum radius (smaller than the former one) and so on. It 
describes therefore an oscillating process where the energy of the explosion decreases each 
cycle, due to the emission of waves during the explosion. 

These numerical simulations have been also computed in bidimensional geometry (cylin- 
drical explosion) and in a one dimensional system (planar explosion). It appears that the 
cavitation process occurs for cylindrical waves but not for planar explosion. An explanation 
from classical results on irrotational flows can be advanced [|r^ : in cylindrical and spherical 
geometries, sound waves always contain both a compressive and a tensile phase, whereas 
a planar compression wave might occur without being accompanied by tension. Then, the 
explosion in two and three spatial dimensions give rise to a negative pressure region where 
the cavitation can take place if the tension is big enough; nothing comparable happens in 
one dimension. 

V. DECOMPOSITION OF THE PROCESS 

The feature of the explosion can be captured by the evolution of the point where the 
density reaches its greatest value. The density pp and the position Rp of this point (called P) 
will be, in particular, investigated. We will indeed focus on large initial density peak, {i.e. 
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(5^ ^ 1); it corresponds typically to an initial pressure peak of more than one hundred bars 
in superfluid helium. In the experiment described above, such pressure cannot be obtained 
by convergent sound waves but only when the primitive bubble collapses [p^]. 

Figure (|]) and figure (|[) show respectively the position and the excess density pp — po of 
the point (P) as a function of time for the explosion listed above (5 = 3 and cr = 1). Both 
graphs exhibit, after a small transient, two dynamical regimes. For small time, (such that 
the density pp ^ 1 ), the behavior is consistent with the following scaling laws: Rp (x \/t 
and Pp — Po oc t~^^'^ (we define u and p such that Rp oc and Pp — po ~ t^^). For different 
values of 6 and a, v varies between 0.46 and 0.51, whereas p varies roughly between 0.5 
and 0.6. At large time (in the figures, for t > 1), we obtain Rp oc Cg ■ t and Pp — Po ~ l/t 
{u = p = 1) , Cg being the local sound velocity. In this 6 and a change, these values 

does not vary significantly (typically between 0.96 and 1). 

The small time regime can be viewed as an explosive- like regime: the shock wave is 
supersonic whereas its density is large (pp ^ 1). For the long time, the system is submitted to 
two effects: one for the bubble dynamics, and one for the shock wave that is now transformed 
in a spherical divergent sound wave (as u = p = 1, pp — pq oc 1/ Rp in fact, as for spherical 
sound wave). The crossover between these two regimes occurs when the shock wave looses 
its supersonic speed (in the simulation shown here, it happens for t ~ 1 unit time). Usually, 
if a bubble is nucleated, it occurs also around the crossover time. 

In fact, the cross-over between these two regimes is still valid, whether a bubble is 
nucleated backward or not. However, as the behavior for Rp does not depend on the initial 
energy (and therefore on 6 and a), the excess density is much more sensible. For long time, 
the shock wave is always behaving as a sound wave (p = 1), although the p = 1/2 property 
for short time is observed only for big enough explosion (see figure ^ for 5 = 8); on the other 
hand, for small initial energy, p is diverging strongly from 1/2. It actually goes continuously 
from p = 1/2 for large energies to p = 1 for small ones (where the shock wave is immediately 
a sound wave). 

As the bubble is nucleated beside the tensile wave, it has its own dynamics. One can see 

10 



clearly on figure |^ b) and c) that the bubble interface moves more slowly that the wave. The 
wave goes away form the center, whereas the bubble grows to its maximal radius through a 
Rayleigh-Plesset like dynamics (see and then collapse under the action of the ambient 
pressure (figure d)). We have defined in fact the radius of the gas bubble Rb as the point 
where the density reaches a typical low value: 

p{Rb) =a- po 

We took a = 0.2 below, although we did check that changing the value of a does not alter 
the results. 

Figure (^ shows the radius of the bubble as function of time (for 6 = 2 and a = 1) and 
illustrates that: 

-the bubble is nucleated at time t*. 

-the growth and the collapse of the bubble are not symmetric. The collapse stage is slightly 
longer than that of the growth. We define tc as the time when the collapse occurs {Rf, = 0). 
-the evolution of the radius for short time satisfies for the growth the relation: 



Rb{t) oc (t-t*) 



as well as for the collapse: Rb{t) oc {tc — t)'^ . 

Last, we computed the ratio between the energy of the secondary peak versus the energy 
of the initial one. We found, depending strongly on the initial shape of the explosion and also 
on the initial energy Eq, that this ratio varies between 0.1 to 0.25. Therefore, a large part 
of the energy has actually been emitted by the explosion and only a small part transforms 
in the bubble. In addition, these ratios has been found to be in reasonable agreement with 



the ratios obtained between the rebound bubble and the initial one in the experiments [IT? 



VI. ANALYSIS OF THE EXPONENTS. 



We will in this section follow the general ideas of explosion |T^. The exponent 2/5 
found for the bubble radius is famous in explosion and comes naturally from a dimensional 
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analysis. Actually, the bubble dynamics involves only a small set of variable: the radius 
r, the time t, the energy Eh of the bubble and the density p^. This is an important point: 
the shape of the bubble interface is only determined by pc as it is the spherical solution of 
SNLS which separates the gas (p = 0) and the liquid (p = |pc) phases. Therefore, only one 
dimensionless variable can be obtained from this set: 

Ef_ 

Pcr^ 

consequently we have obtained the well-known dependence of the radius of the bubble 
with time: 

where /? is a constant (determined in fact by the shape of the interface). 

A contrario, before the bubble is nucleated, the shock wave obeys different scaling and 
such a dimensionless variable is not unique. Indeed, the density of the shock waves also 
varies now. But, as figure shows, for huge explosion (here 5 = 8), a self-similar shape 
in the density profiles can be observed. We will therefore limit our investigation to big 
explosions (5^1) and to time such that p ^ 1 everywhere inside the shock region; in 
addition, the following remarks are needed: 

-first of all, a self-similar solution cannot conserve both energy and mass, (this follows 
from the fact that mass is an integral of {ipl"^, whereas energy is an integral of IV'I^)- 

-consequently, the limit of infinitely thin explosion is either finite energy with zero mass 
or finite mass and infinite energy. 

- due to the condition at infinity, (p = po), a flux of mass is continuously entering the 
self-similar region as the time goes on, whereas such mass density has zero energy density. 

Balancing terms in Bernoulli equation (H) motivated us to look at solutions such that 
(the cubic term is neglected versus the quintic one in SNLS): 

p(r-t) = l=h\0 <P{r,t)=giO with ^ = (6) 
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In this approach, ^ is the self similar variable. It has the same structure as the one 
elaborated for the heat equation. In both case, it comes from the balance between the time 
derivative and the laplacian term. Therefore, the position of the peak of density is defined 
hj ^ = a, where a is determined by the initial condition and is time independent. Equation 
(^) is valid everywhere except that because of the assumption (p ^ 1), it is relevant and 
consistent only near the shock wave region. Also, in our regime of interest, we have a ^ 1. 

It follows from this analysis that: 

7-. / N r 1 h'^(a) 
Rp{t) = ayt and pp — pa = — -j=- 

yt 

which are in good agreement with the numerics. Notice that if the quintic term were 
negligible versus the cubic one in SNLS, the scaling for the density Pp — po would have been 
in 1/t. It explains why when the energy of the explosion is lowered, the value of u goes 
continuously from 1/2 to 1. 

Finally, with these scaling laws, the mass inside the shock sphere evolves linearly in time, 
whereas the available flow of particles goes like the square root of time. Then with respect to 
mass conservation, these scaling laws cannot apply as time goes to zero, when the entering 
flow of particles cannot balance the growth of mass of the solution. This gives a critical 
time under which the self-similar solution is not valid. Numerically, it is difficult to see if 
the transient that we observe initially is due to this constraint or to numerical ones. 

The radial velocity v is obtained through: 

r 

and the Bernoulh equation (||) reads: 

Another equation is needed for describing the dynamics; the natural one would be the 
mass conservation. In fact, the integral equation for the conservation of energy in a sphere 
of constant radius in the self similar variable = constant) is more useful: 
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V ■ ((Vv^)' + Ipv' - ^v^AVp + p') + ^VpV(v/p)A0 = ^R'it)eip,v) (8) 
e being the density of energy: 

Equation (H) reflects the fact that the increase of energy inside the sphere of radius ^ is 
balanced by the enthalpy flux. 

An exact solution of the system can be found by neglecting the quantum pressure: 

^ and p = \ - X ^ 



4Vi \ 2 AVi 

corresponding to 5^(0 = ^^(0 = ~ ^ ~ (9) 

8 V 2 4 

These solutions also satisfy exactly the continuity equation (^). Identifying the energy 
of such a solution with the initial energy Eo gives: 

E, = = l/a« 

° 15 64V 2 

and therefore: 

a oc ^/o^S 

Figure a) shows the phase as function of ^ at different times for 6 = 8 and ct = 1; it 
offers a good agreement with the parabolic analytical solution of equation (^). Particularly, 
the minimum of the curve (at ^ = 0) does not change its value as time goes on, as predicted 
by the self similar solution. In addition, figure (|^ b)) shows for different 6 (with a = 1), 
(t>{C,/a)/a^; the solution (|^) suggests that these curves should coincide. The numerical results 
are less accurate than the one for a given 6, although they are still reasonable. 

On the other hand, the density profile provided by the solution (P) does not agree well 
with the numerics (see figure ^ for 6 = 8): although the linear behavior for p is acceptable 
in an intermediate range (1 ^ ^ ^ a), there are strong differences for ^ ~ and ^ ~ a. For 
^ < 1, the quantum pressure becomes dominant in the self similar set of equations (|^ and 
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(P). For ^ ~ a the solution is in fact maximal and also has to match with the condition at 
infinity. 

As the numerical simulations show a good agreement with the self similar solution for 
the phase, we will assume the solution (|^) to be exact for cj) and therefore we will try to 
find a more accurate solution for p, particularly in restoring the quantum pressure. After 
rescaling h and ^ for convenience {h{i) = 1^/(31/4^/2)), the problem reads: 

nu) + 2^ = f{u)-u'f{u) (10) 

The boundary of the self-similar region being then defined by u = a' = 3^^^a/2. Even if 
the equation is relevant only for u < a', it is interesting to study it for all values of u. Using 
Mathematica |1^, it is straightforward to analyze the class of numerical solutions of (p!OD 
with the physical initial conditions /(O) = A and /'(O) = 0. Depending on the value of A 
compared to a critical value Ac = 0.98473..., we obtain the following behaviors: 

-for A < Ac the solution oscillates and goes to zero at infinity. 

-for A > Ac the solution includes a singularity in the real plane. 

Figure (|^) shows the solution /(m)^ as a function of u for two values of A: one slightly 
smaller than Ac, the other slightly larger. Ac has been evaluated through the shooting 
method. The solution for A = Ac will behave linearly in m as m — > 00, in agreement with 
the solution (H). 

The solutions for A < Ac appear to be in good agreement with the numerical solution 
of SNLS (see again figure H), and we will therefore assume that it describes the self similar 
solution. The exact value of A will be a function of a and consequently of Eq. In addition, 
the matching between this solution around u ^ a' with the sound wave solution of SNLS 
(at density po) inform us about the wave number of the oscillations emitted by the solution. 
The oscillating solution at large u for A < Ac can be obtained easily by neglecting when 
u ^ oo] then the equation: 

u 
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has an exact solution in terms of the Bessel functions of the first kind Jr,(-)- It reads: 



f{u) = 



Ci and C2 being two constants. 

For large x, the behavior of J-|-i/4(x) is simple: 

cos{x + <^) 
J±i/a{x) 

where <y9 is a constant. 

Therefore we obtain the following behavior for f{u) for large u (more precisely for m 3> 1): 

/(") « 

So that, in the neighborhood u ~ a', we have: 

, cos{a' ■ Sx + ip') 
+ '^^) ^ 

This corresponds to oscillation for the density p, around r = a\/t, with wave number 

k{ty. 

This simple matching between the self similar solution and the density at infinity through 
excitation waves gives in fact the wave number of these radiations as function of time. 

VII. RENORMALISATION GROUP CALCULATIONS. 

Through renormalisation group (RG) calculations one can describe precisely the 
behavior of f{u) at both edges of the shock wave u = and u = a'. 

for ^ — s> 0: introducing an arbitrary small parameter e, through rj = ^, the following 
equation for m{ri) = f{u) is obtained: 
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m"{7]) + = e^m'^iT]) - e^r]^m{r]) (11) 



Then m{ri) is expanded as a series in e^: 

m{r]) = mo{r]) + e'^mi{r]) + e^m2{ri) + .. 
The system of equation for mo and mi reads: 

rriQ + 2^ = 



ri 



m'l 



m'i + 2— = mo^ 
?7 



The boundary conditions for m come from the constraint m'(0) = (no flux at = 0) 
and from the asymptotic behavior at large r], which will act as a matching condition for the 
remaining constant of integration. 

Using m'(0) = 0, one obtains: 

2 2 

mil]) = A(0) + ^A(O)^ 

yl(0) being the constant of integration at t = 0[|. 

This formula loses its validity for er/ > 1 because the perturbative expansion of m is 
meaningless, (e^mi is no longer a small correction to mo in this regime). 

The renormalisation group (RG) method consists in considering the integration constant 
as being slightly dependent on rj through the change: 

A(0) = Z(0, ii)A{mu) with Z(0, /i) = 1 + ^ an(yu)e^" 



^formally, the other solution of this second order differential equation m = B/rj should be kept for 
the formal solution of the problem because the boundary condition m'(0) = has to be imposed 
only for the final results; however, as we checked that this term has zero contribution at the order 
of the analysis, it has been forgotten for convenience. 
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Z is called the multiplicative renormalization constant and a„ is chosen such that the solution 
m{r\) has the same structure for all 

This imphes that: 

Obviously m{ri) should not depend on /j,. Then, the so-called RG equation d^m{rj) = gives: 

dA 1 



d/j, 3 ^ 



whose solutions read: 

Ail,) 



6 V^' 



Substituting A{fi) into the expression for m{7]) and putting fj. = rj gives the following 
formula for the initial function /: 

/ 6 V^' 

f{u)- 



Xo is a constant of integration to be determined. Notably, we have obtained the kind of 
divergence that occurs for A > Ac- 

For ^ — > a: the same treatment works in the limit a' ^ oo by writing: 



f{u) = Va' + — 



(1/a' has the same feature that e had formerly and the l/(a')^ factor has been introduced 
in order to neglect the non linear terms for the first order correction). 
Eventually, it gives the following behavior (Aq being a constant) : 

J [u) = V a' — . e 8i' cosh(a [u — a) ) 

[a) 4 

Aq and Xo are yet to be determined. This is usually accomplished through a matching 

condition written in an intermediate region between these two regimes. In our special 
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problem, an easy way to perform this would be to match both solution in the region 1 ^ m ^ 
a' where we know that the solution reads f{u) ~ ^/u. In fact, the first order RG calculations 
is not accurate enough to realize such a matching; but the goal of these calculations was 
actually to have a reasonable idea of the shape of the solution near both edges. The whole 
matching process would need a more complete analysis of the RG theory. 



VIII. CONCLUSION-ACKNOWLEDGMENTS 

Finally, as the shock emits waves, the energy of the self similar solution cannot be 
constant in time. Also, the cubic term has been neglected for finding the above self similar 
solutions. This should really be incorporated into a pertubative analysis by considering a 
(or equivalently Eq) as being slowly dependent on the time. Then, through a solubility 



equation, an evolution equation should be found for a [|20|. This detailed analysis, which 
would provide for example the amplitude of the emitted waves, will be the subject of a 
further work. 

It is a pleasure for me to thank Yves Pomeau, Leo Kadanoff, Shankar Venkataramani 
and Alberto Verga for helpful discussions and for their interest on this paper. This work 
has been supported in part by the ONR grant: N00014-96-1-0127 and the MRSEC with the 
National Science Foundation DMR grant: 9400379. 
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FIGURES 
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FIG. 1. Respectively, the pressure as function of the density and the sound velocity as function 
of the pressure for the model SNLS. These graphs show a consistent liquid-gas transition, where 
the spinodal decomposition point is located at p = pc- 



22 



5 




10 20 30 40 10 20 30 40 

radius r i \ radius r 



FIG. 2. Density profile for (5 = 3 and cr = 1 at various time a) for small time, t=0.02, t=0.06 
and t=0.12 unit time; the explosion gives rise to a shock wave, b) for t=0.18 and t=0.6; a bubble 
appears beside the tensile part of the shock wave; oscillating perturbations are emitted by the 
shock wave, c) for t=2.9, t=7.9 and 13.4; the bubble grows until a maximal radius and begins to 
retract at t=13.4; the shock wave has changed in a spherical sound wave and is no more moving 
with the bubble interface, d) t=18.4, t=23.4 and t=24.4 the bubble collapses, giving rise to a 
secondary explosion. A large part of the energy has been radiated through the spherical waves. 
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FIG. 3. Position, as a function of time, of the maximum of the density. It exhibits two typical 
behavior, one for small time, when the wave is supersonic, the other one later, where the shock 
wave has transformed in a sound wave (the position evolves then linearly in time). For the small 
time, the position of the shock wave is consistent with a law Rp oc t^^"^ . 
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FIG. 4. The deviation from average density as function of time. After a short transient, this 
quantity evolves as a power law in time {pp — 1 cx t~'^). The two lines above the graphs correspond 
to /X = 0.5 for small time and = 1 later. 
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time t 

FIG. 5. The over-density as function of time for a huge explosion 5 = 8 and a = 1. The hne 



indicates the slope t The small time behavior coincides with a /x = 1/2 regime. 
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FIG. 6. Radius of the bubble induced by an explosion as function of time {5 = 2 and cr = 1). 
The growth and the collapse are not symmetric. Both the growth and the collapse obey the 
Sedov-Taylor law for the radius as function of time. 
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FIG. 7. a) phase ^ as a function of ^ the self similar variable, for different time (from t = 0.01 
to t = 0.2; 6 = 8). Except the two first curves, which correspond to the lowest time, where the 
dynamics is still in the transient, the different curves almost coincide for ^ ~ and have a parabolic 
shape, as predicted by the simple self similar approach, b) for different 6, (f)/6'^ as function of 
according to the self similar analysis, these curves should coincide. 
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FIG. 8. density profile of the solution of SNLS for 5 = 8 and a = 1 for different time smaller 
than 0.05 unit time. One can observe a self similar shape in the shock wave as time goes on. Even 
if the density obeys a linear profile for intermediate values of r, it is not valid anymore for r — > 0. 
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